Resid0 = NaN(size(Y_reg,1),1);
Resid_ev0 = NaN(size(Y_reg,1),1);
%[iN jN] = find(isnan([Y_reg X_reg]));
%iN  = find(any(isnan([Y_reg X_reg]')));
iN  =  any(isnan([Y_reg X_reg]')); 
Y_reg(iN,:) = [];
X_reg(iN,:) = [];

%% eliminate columns that are all zero
jj_col = all(X_reg==0);
X_reg(:,jj_col) = [];
if shrink
  X_d(:,jj_col) = [];
end

jj_col_l = find(jj_col);
if any(jj_col_l <= ev_loc)
    ev_loc0 = ev_loc-length(find(jj_col_l <= ev_loc));
else
    ev_loc0 = ev_loc;
end

jj_col_l = jj_col_l((jj_col_l<length(reg_names_reg)));
reg_names_reg(jj_col_l) = [];
 

nobs = size(Y_reg,1);
if ( rank(X_reg) == size(X_reg,2) ) & (rank(X_reg) > 0 )%% nobs > size(X_reg,2)
  
  iXX_d = inv(X_reg' * X_reg + X_d'* X_d);    
  beta = iXX_d*(X_reg' * Y_reg);
  %  beta0 = (X_reg' * X_reg) \ (X_reg' * Y_reg);    disp([beta beta0]);  pause;
  resid = Y_reg - X_reg * beta;
  s2 = (resid' * resid)/nobs;
  var_ols = s2*iXX_d;
  
  if ispanel
    %% var/cov matrix of residuals
    Resid0(~iN) = resid;
    Resid0 = reshape(Resid0,[tT nN]);
    
    if sandwich
      %% available obs for each forecaster
      nNobs = sum(~isnan(Resid0),1)';
      mnNobs = min(repmat(nNobs,[1 nN]),repmat(nNobs,[1 nN])')  ;

      Resid00 = Resid0; 
      Resid00(isnan(Resid00)) = 0;
      VarR =  (Resid00'*Resid00)./mnNobs;
      %[nanvar(Resid0,1)'+nanmean(Resid0,1).^2' diag(VarR)]
      

      if sandwich == 1
          %% diagonal
          Sigm = diag(diag(VarR));
          %% check whether works w OLS % Sigm = s2*eye(nN);
    
      elseif sandwich == 2
              %% principal component -- equal loadings
              aggRes = nanmean(Resid0,2);
              %% var common
              varagg = nanvar(aggRes,1)+nanmean(aggRes).^2;
              Resid00 = Resid0 - repmat(aggRes,[ 1 nN]);
              %% variance idiosyncratic
              VarR = nanvar(Resid00,1)'+nanmean(Resid00,1).^2';
              Sigm = diag(VarR)+varagg*ones(nN,nN);
              
          elseif sandwich == 3
              %% principal component
              [Vv,Dv] = eig(VarR);
              Sigm = Dv(end,end)*Vv(:,end)*Vv(:,end)';
              %% put back originial diagonal
              Sigm(find(eye(nN))) = diag(VarR);
          end
          
          Sigm = kron(Sigm,eye(tT));
          Sigm(:,iN) = [];
          Sigm(iN,:) = [];
          SigmX = Sigm*X_reg;

          XpSigmX = X_reg'*SigmX;
            iXX = inv(X_reg' * X_reg);    
      var_beta = iXX*XpSigmX*iXX; 
      
      if Gls
          iSigmX = inv(Sigm)*X_reg;
          XpiSigmX = X_reg'*iSigmX;
          XpiSigmY = iSigmX'*Y_reg;
          beta = XpiSigmX \ XpiSigmY;
      end
      
      
    elseif ~sandwich
      var_beta = var_ols;      
    end    
    
  elseif ~ispanel
    var_beta = var_ols;
  end
  
  resid_ev0 = Y_reg - X_reg(:,ev_loc0+1:end) * beta(ev_loc0+1:end);
  stdev = diag(var_beta).^.5;

  fprintf(fileID0,' nobs: %2.0f; s2: %2.4g, s2/varY: %2.4g\n ',[nobs,s2,s2/var(Y_reg)]);
  for jj_r = 1:length(reg_names_reg)
    fprintf(fileID0,' %2.19s: ',[reg_names_reg{jj_r}]);
    fprintf(fileID0,' %2.4g (%2.4g) [%2.4f]; ',[beta(jj_r) stdev(jj_r)  beta(jj_r)/stdev(jj_r)]);
  end
  
  info = [s2,s2/var(Y_reg),nobs];

  if ~isempty(iN)
    Resid_ev0(~iN) = resid_ev0;
  else
    Resid_ev0  = resid_ev0;
  end

else %% not enough obs
  fprintf(fileID0,' nobs: %2.0f;  not enough obs \n ',[nobs]);
  
  beta = NaN(size(X_reg,2),1);
  stdev = NaN(size(X_reg,2),1);
  info = [NaN,NaN,nobs];

end

